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In this work we study numerically and analytically the interaction of two qubits in a one¬ 
dimensional waveguide, as mediated by the photons that propagate through the guide. We develop 
strategies to assert the Markovianity of the problem, the effective qubit-qubit interactions and their 
individual and collective spontaneous emission. We prove the existence of collective Lamb-shifts 
that affect the qubit-qubit interactions and the dependency of coherent and incoherent interactions 
on the qubit separation. We also develop the scattering theory associated to these models and prove 
single photon spectroscopy does probes the renormalized resonances of the single- and multi-qubit 
models, in sharp contrast with earlier toy models where individual and collective Lamb shifts cancel. 


I. INTRODUCTION 

Technological progress in the fields of superconducting 
quantum circuits [T] and nanophotonics [2] are opening 
the door to new quantum technologies based on propa¬ 
gating photons and few-level systems, such as few-photon 
transistors US], non-classical states of the radiation field 
[SHH or new photodetector schemes [30. Of particular 
fundamental interest are also proposals to use these low¬ 
dimensional photonic systems to mediate interactions be¬ 
tween qubits [SHn], for which there are already initial im¬ 
plementations in the superconducting world HI ESI, and 
which could potentially have applications in the quantum 
information world, as well as in quantum simulation and 
the study of quantum phase transitions. 

At the same time that this technology evolves, we 
need to further develop the theoretical tools that study 
the light-matter and light-mediated matter interactions. 
While there exist approximation methods for some of 
the optical technologies mentioned above, based on 
single-photon single-qubit effective boundary conditions 
[HEl], input-output theory [Him or scattering theory 
mm, the degree of control of some methods or their 
computational generality for any number of photons is 
an open problem. 

Matrix Product State (MPS) methods are one novel 
tool that solves some of the limitations of earlier meth¬ 
ods regarding photon numbers and interaction strengths 
to describe the full light-matter wavefunction and its time 
evolution. These methods rely on a discretization of the 
photonic degrees of freedom either in frequency or in po¬ 
sition space, and seem to provide qualitatively and even 
quantitatively accurate results for problems of few pho¬ 
ton scattering [TOMT] . Two questions open: how can 
we validate these models and discretizations and how do 
they compare to the usual techniques, i.e. master equa¬ 
tions and Markovian approximations. 

This work dwells into these problems by studying a 
very precise feature of the numerical approximations 
mentioned above: the discretization of the photonic 
degrees of freedom. We study the effective interac¬ 
tion of two qubits in the weak- and strong-coupling 
regime, both within the Rotating Wave approximation 


and with Matrix Product States, for a variety of un¬ 
derlying microscopic models of the photonic degrees of 
freedom. Combining theoretical estimates based on the 
Wigner-Weisskopf theory with exact analytical studies of 
the single-excitation or Rotating Wave Approximation 
(RWA) limit we prove that the Markovian approxima¬ 
tion is only well justified for qubits inside their mutual 
light-cones and develop accurate numerical procedures 
to extract the effective parameters of the qubit-qubit in¬ 
teraction and collective decay. These parameters agree 
qualitatively with earlier works based on the resonant- 
dipole approximation 0. However, we also find collec¬ 
tive effects (Lamb-shifts) which strongly depend on the 
cut-offs and the microscopic details of the interaction, 
and which renormalize the qubit frequencies, their inter¬ 
actions and the effective individual and collective dissipa¬ 
tion. The renormalization effects can be of measurable 
magnitude, increasing with the qubit-qubit separation, 
and we suggest how recent experiments with supercon¬ 
ducting circuits could probe those features. 

There are two reasons why one could doubt of these 
non-universal renormalization effect: such effects have 
not been found in theoretical studies of single-photon 
scattering [12 EH and we have derived them using RWA 
equations. We dispel the first criticism studing the 
Lippmann-Schwinger scattering formalism and proving 
that its resonances are determined by the same renormal¬ 
ized parameters of the theoretical model. We also study 
existing single-photon scattering literature and conclude 
that the ’’toy”-models used to represent the photons are 
often too simple and have zero Lamb shifts. Finally, we 
report on quasi-exact numerical simulations of the full 
spin-boson model including counter-rotating terms, for 
exactly the same parameters of the RWA simulations. 
These MPS simulations confirm the RWA results for 
coupling strengths below five percent of the qubit gap, 
g/A. < 5%, which is the current operating regime for 
circuit-QED in open transmission lines [312 (HUH. 

The structure of this manuscript is as follows. In 
Sect. [^ we introduce a simplified light-matter interac¬ 
tion model for two qubits that talk to a one-dimensional 
photonic medium. We derive the evolution equations for 
the qubits in the single-excitation limit and obtain the 


2 


Markovian approximation and resulting master equation 
with care. This derivation pays great attention to the 
infrarred and ultraviolet cut-offs of the problem, leading 
us to conclude that individual and collective Lamb-shift 
influence the qubit frequencies, and the spontaneous and 
collective emission rates of the qubits beyond earlier pre¬ 
dictions. In Sect. |IV| we relate the parameters of the 
Markovian model with the Lippmann-Schiwnger scatter¬ 
ing states and dicuss why earlier scattering studies did 
not show the renormalization effects. While the theoreti¬ 
cal results so far were based on a few solvable models, 
in Sects. |VB| and |V C| we return to physically realis¬ 
tic models of waveguides and of photonic crystals. We 
confirm the role of the cut-off frequency in the single¬ 
qubit renormalization and in the two-qubit interactions 
in both cases. In Sect. fVll we introduce numerical simu¬ 
lations based on Matrix Product States and confirm that 
they agree with the theoretical predictions in the strong¬ 
coupling limit. Finally in Sect. m we summarize the 
main consequences of this study and discuss possible ex¬ 
perimental studies of these non-universal effects. 

II. THEORETICAL MODELS 


Using the fact that momenta are equispaced, and that the 
frequency spacing dw = uj'{k)dk = dk/p(uj) is inversely 
proportional to the density of states, p(a;), we may write 

J(a;) = 2 X 27r / g‘^ujk5{uj - uJk)piyJk)duJk (5) 
Jo 

= ing'^ojp^oj). 

If we assume the RWA and study the spontaneous emis¬ 
sion of a two-level system interacting with this waveg¬ 
uide, we will obtain that the spontaneous emission rate 
is 7/2 = J(a;)/2 = g'^ujp{uj), independent of the micro¬ 
scopic details of the discretization we used. 


B. Microscopic models 

We will typically consider two types of problems: gap¬ 
less photons and a photonic crystal with a finite band¬ 
width. The first case corresponds to a chain of coupled 
oscillators and it models the propagation of optical pho¬ 
tons or microwave photons in a waveguide |I9j . The dis¬ 
persion relation has the form 


A. Spin-boson model 


We wish to study the interaction of one or more two- 
level systems with a photonic medium. For that we rely 
on the spin-boson model, which is a convenient descrip¬ 
tion for the type of implementations where this problem 
has some interest —namely superconducting circuits and 
nanophotonics. The model reads 

i k k,i 


The light-matter coupling typically adopts the form 


9ik 





( 2 ) 


where uJk is the dispersion relation, Xi is the position of 
the f-th qubit in real space, L is the physical length of the 
one-dimensional medium and 5 is a parameter measuring 
the coupling strength. 

As discussed below, the model § adopts different dis¬ 
persion relations and coupling strengths depending on 
the microscopic description of the underlying photonic 
medium. For computational reasons, in any of these de¬ 
scriptions we will have to impose both a finite medium 
size, L, and a spatial discretization Aa; = L/N. This 
leads to a discretization of the momenta 

A: G — X {0,±I,±2... ± y}. (3) 

In either case we may introduce the spectral function 


J{uj) = 27r'^\gik\'^ S{uj - Wk). (4) 

k 


uJk = uJcVJ - 2 cos(fcAa;), gsk = ( 6 ) 

with a cut-off Wc = 1 /Ax that depends on the discretiza¬ 
tion. With this choice of units, the dispersion relation 
becomes linear at low momenta, ujk — v\k\, with unit 
speed of light, v = 1. 

The other possible problem that we will study is a 
chain of coupled cavities. In this case the dispersion re¬ 
lation has a finite bandwidth 


Uk = ojo- J cos(fcAx), gsk = g\j (7) 

and the discretization step plays no role, so that we can 
take Ax := 1. 

Finally, in addition to these models with a finite num¬ 
ber of modes, we will consider a continuum limit in fre¬ 
quency space that has L —>■ 00 and introduces the cut-off 
in the coupling strength, not on the dispersion relation 

u^k = 1 * 1 , 9sk = ( 8 ) 

This exponential cut-off is usually employed in studies of 
the spin-boson model [H] because it allows efficient com¬ 
putation of memory functions and analytical predictions. 

In addition to these realistic models, Sect. |IV C| an¬ 
alyzes other models that have become standard in the 
literature, such as the linearized model by Shen and Fan 
m [T5] or coupled cavities models [23l |24j, where the 
coupling strength is approximated as frequency indepen¬ 
dent, leading to J{uj) approximately constant. 
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C. RWA and single-excitation equations 


In the low-energy and small-coupling limit, in which 
g is small compared with the qubit and photon frequen¬ 
cies m, we can simplify the spin-boson model adopting 
the so called rotating-wave approximation, in which the 
counterrotating terms -I- H.c. are neglected. In the 

RWA limit the number of excitations in the system is 
conserved and we can write down variational wavefunc- 
tions in which either one of the qubits is excited, or a 
photon is propagated. We write this as 




C,cr+ + 

. i k . 




(9) 


decays much faster than the time it takes for the qubit 
to evolve. 

If we have a single qubit, we will be concerned with the 
evolution of the single memory function Kss{t). We as¬ 
sume a-priori, that c(t) evolves as exp{—iA't)w{t)^ with 
some renormalized frequency A' that is close to the orig¬ 
inal qubit frequency A, and a slow function w(t): 

bss(t) =-i f Kss(t - T)e~^^ ^w{T)dT (15) 

Jo 

f Kss{t — 

J —OO 


This leads to a memoryless equation for one qubit 


where jn) = ®i \gi) | 0 fe) is the vacuum in the qubit 


and photon spaces. 

The evolution equations for the state read 


'^dtCs — ^ ^ 9'i'k'^k j 

k 

(10) 

idtipk = uJk'ipk + Y 9*skCs- 

S 


We integrate formally the photons and plug this expres¬ 
sion into the qubits 

T ds (^) H" ^ ^ ^ss' (j') ■ 

s' 

(11) 

Each qubit has an input signal that either comes from the 
photonic field, or from a memory function that concerns 
the dynamics of all of the qubits in the past 

as{t) =Y9ske *‘^''‘7fc(0), 

k 

(12) 

bss' (t) = -M Kss' {t - r)cs' (r)dT 

Jo 

(13) 


with the kernel 


idtCg = (a), - ly) Cs, (16) 

replacing the memory function with the constants 

ry 

6g-i^= 2A(r;0)e*^>dT. (17) 


Using the identity 

pOO 

Jo 


tS{uj - A') - *PV — 


A' -w 


(18) 

we obtain the correction to the qubit frequency and the 
decay rate as solutions of the self-consistency equation 

= = (19) 

= JiK) (20) 


In Fig. we summarize the outcome of these compu¬ 
tations for the exponential cut-off (|^, for which we can 
compute the kernel exactly 


Kgg^it) 





^—lk{Xs —Xg/ ) 


K{t;dss') + K{t]-dss')- 


c.c. 


dwfc 


(14) 


-f = g‘^/S!e 

r 2 1 r 




Ei(A7a;c 


( 21 ) 


For a linear dispersion relation, Wk = v\k\, the previous 
kernel adopts the form K{t + dss')+K(t — dss'), implying 
that the separation between qubits induces a delay in the 
mutual interaction. 


III. MASTER EQUATIONS 
A. Single-qubit master equation limit 


where Ei(x) is the exponential integral, a negligible term 
when Wc —> oo. Note how the renormalization of the 
qubit frequency 5 may become significantly larger than 
the decay rate itself, diverging as the cut-off increases. 
However, even if this non-universal effect is compara¬ 
tively large, the ratio <5/A is still small and Fig.[^ shows 
that we can solve Eq. (19) in a simpler, approximate way, 
using A instead of A' in the integral, and replacing the 
resulting correction A' in the expression for 7 


Equation (10) can be manipulated to study the dy¬ 
namics of the qubit and the field under further approxi¬ 
mations. The most common one is the Markov approxi¬ 
mation in which we assume that the memory of the field 


(A^^_ydw, (22) 

7. = J{K)- 
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FIG. 1. Parameters for the evolution of the two-level system 
coupled to the ID transmission line as a function of the cut¬ 
off, computed with Eq. (171 for model (Hf. In (a) we fix 
A = 1 and try different cut-offs. In (b) we fix the cut-off 
ujc = 10A and change A'. It is evident from this plot that 
since |A — A'| — <C 1, there is no big difference between 
using A and A' in Eq. |T^. 


B. Two-qubit master equation limit 

The same procedure can be repeated for two or more 
qubits. In this case the outcome should read 


idtCs = 








with the additional constants 


. ^ss' I 

9ss' - = J Kss'[T)e 


iA' / r 


dr. 


(24) 


There are several caveats on Eq. (231, though. The first 


one is that this equation now is only valid for times 
t > Tss' = Idss'l/t'; larger than the time it takes a photon 
emitted by qubit s to arrive to qubit s'. Before that time, 
as we will see below, the qubits evolve almost indepen¬ 
dent from each other, in an approximate causality. 

There are also some subtleties when we try to estimate 
the parameters gss' and 7 ss'. Firstly, Eq. (24) neglects 
the effect of the collective coupling gss' in the renormal¬ 
ization of the single-qubit frequencies Ag. Secondly, gss' 
is related to the principal value part of the spectral func¬ 
tion. We cannot compute this in closed form because the 
RWA precludes the use of the Kamers-Kronning relation, 
but assuming that non-RWA terms are negligible we can 
approximate gss' by the expression that results from the 


full model (See Suppl. Matt, in [^). Based on these 
considerations, an approximate set of constants for this 
problem that is often used in the literature, combines Eq. 
(19) with the resonant dipole approximation 


7ss' — 7s COs(/Csdss') 7 dss' 

7^ 

Qss' - y Sin(fcsldss'l)) 


where kg > 0 is the positive momentum associated to 
A(,. An important prediction of this formula is that, as 
long as it is valid, dipole-dipole interactions can he sup¬ 
pressed by placing quantum emitters at distances such 
that /csids.s'l = WTT. This condition leads to the imple¬ 
mentation of homogeneous superadiant models [3I1S]. 

When the qubit have identical frequencies, Ai = A 2 = 
A, we can develop a more accurate approach to the 
Markovian limit that takes into account self-consistent 
renormalization. In this case the exact problem decou¬ 
ples into two equations, with variables c± = ;^(ci ± C 2 ) 


idc± = Ac± — i K±{t — r)c±(r)dT ~ (A± — i^±)c±. 
Jo 

(26) 

with K±{t) = 2K{t\ 0) ± [K{t] d) -|- K{t; —d)] and d the 


distance between the two qubits. From these equations 
we recover 


A' = i(A+ + A_), 
gi 2 = ^(A+-A_), 


7 = 

712 = - 1 -) 


(27) 


in a more correct than the derivation above. Let us re¬ 
mark the need to compute A± separately, using either 
(19) or ( [ 2 ^ . This gives rise to collective renormalization 
effects that depend on 512 and these effects influence the 
values of 7 and 712 , in contrast with Eq. (25), where the 
spontaneous emission parameters only depended on A'. 

We can particularize earlier expressions for the model 
([^. Keeping the lowest order terms in the coupling, g, 
and dropping the corrections from using Eq. (27), the 
Markovian interaction reads, 


A'= A — 5^ —[wc + e ^/‘^“ifi(A/wc)], 

TT 

g^ “Jojc 
27r 1 -|- d^ujc 

+ -(sin(dA) + /(dA)), 

7 = -b ©(/wc) 

712 = cos(dA) -b ©(g'^Wc), 


(28) 

(29) 


(30) 

(31) 


where we introduce f{x) = —Re -b A/wc)]}. 

This result agrees qualitatively with the approximation 
in Eq. (25) in the limit dA ^ 1, where interaction be¬ 
comes approximately periodic. The qubit interaction, 
however, diverges at short distances and becomes of order 
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FIG. 2 . (a) Parameters for the evolution of the two qubits 
coupled to the ID transmission line as a function of the dis¬ 
tance between qubits, measured in units of Ao, the qubit 
wavelength. We plot 7/2g^ (dash-dot), 712/2^^ (dashed) and 
gi2lg^ (solid), computed using \ 27 \ , using a self-consistent 
formula for A±. (b-d) Errors in the parameters 5112, 712 and 7 
made by using no Lamb-shift renormalization (dashed) or by 
using the simplified renormalization scheme ( 22 1 (solid). Note 
that the error is amplified with respect to the single-qubit case 
because already 7 ~ 271. In all simulations ujc = 10 A and 
g = 0 . 04 A. 


6 = —g^ujc- This happens at a distance that is inversely 
proportional to the cut-off, which is the quantity that de¬ 
termines our spatial resolution. Moreover, as we will see 
below, the renormalization of the frequencies introduces 
large discrepancies with ( 251 . 

In Fig. we plot the constants that result from an 
exact self-consistent evaluation of the Markovian param¬ 
eters through a direct evaluation of the kernel integral for 
the model (|^. The inclusion of the self-consistent correc¬ 
tions does not change the qualitative shape of gi2 and 712, 
but because A+ is not identical to A_, there are beatings 
that appear in 7 and which amplify with the distance and 
the cut-off. In order to investigate this further, Fig. we 
have taken the difference between the full, self-consistent 
solution in Fig. [^, and a solution without Lamb shift 
corrections (dashed) and where these corrections are in¬ 
troduced only in 7 (solid), as in Eq. ( 22 ). These differ¬ 
ences are shown in Figs. -d for gi2, 7 and 712. Unlike 
for a single qubit, there are sizeable discrepancies even 
in the coherent part of the evolution gi2. Note also how 

7 differs enormously from the trivial, distance indepen¬ 
dent prediction, and presents those beatings mentioned 
before. Finally, 712 also differs from the unrenormalized 
case, but it seems that already Eq. provides a good 
enough approximation. 


IV. LIPPMANN-SCHWINGER SCATTERING 
THEORY 


So far we have derived effective models for the nonequi¬ 
librium dynamics of the qubits themselves, tracing out 
the photons. We will now prove that the resonances of 
this effective master equation may be probed using single¬ 
photon spectroscopy to excite the internal transitions of 
the effective spin-spin interaction and accessing 7, 712 
and gi 2 through the changes in state of the photons that 
interact with the qubits themselves. 

Given the constraints of our model, which is developed 
in frequency space, we cannot resort to the direct compu¬ 
tation of scattering matrices [TillTS] . Instead we will use 
the resolvent method, introducing the basic idea in Sect. 


IV A and applying it to single- and two-qubit processes 
in subsequent texts. 


A. Scattering states 

We study the propagation of waves when they face 
an impurity or potential represented by a local interac¬ 
tion V following the operator formalism that was set up 
by Lipmman and Schwinger and is nicely summarized in 
Ref. [ 26 ) . The original problem admitted a continuum of 
solutions labeled by some momentum k and energy Ek 

{Ho - Ufe) |$fe) = 0 . ( 32 ) 

The new solutions are expected to assimilate to the re¬ 
flected beams out of an incoming wave |<i>fe). Because far 




















































6 


away from the perturbation they must contain such plane 
wave, and at those distances the dynamics is dominated 
by i?o, these solutions must have the same energy E but 
satisfy a different equation 

{H-Ek)\^k) = {Ho + V-Ek)\^k)=Q- (33) 


within the single excitation sector. Let us introduce the 


right-hand-side of Eq. 35 


\4') = vm = 




If7), 


(39) 


For simplicity we are also going to drop the k label unless 
explicitly needed. 

The Lippmann-Schwinger (LS) formalism starts by 
defining the resolvents of the free and interacting prob¬ 
lems 


Go{E) = {E- Ho)-\ G{E) = {E - H)-\ (34) 

We also introduce the quantity E^ = E + ie, modifying 
the scattering problem to 

{E+- = 0. (35) 


where e > 0 will be eventually taken to the limit e —> 
0+. The sign of this small perturbation is relevant to 
distinguish the direct problem ($ represents the incoming 
wave), from the inverse problem ($ is the time reversed 
of the scattered wave). 


From (35) we construct a self-consistent solution 


|4/) = |4>) + Go{E+)V Id/) = J2^Go{E+)vr |$) (36) 


n—0 


If we apply the operator (if+ — iJo) on this equation and 
use the fact that |$) belongs to its kernel, we recover 
(35). The 1$) provides the appropriate incoming bound¬ 


ary condition, while the second term represents all the 
scattered waves. 


The LS equation (36) can be solved self-consistently. 


The first order truncation of the series is the so called 
Born approximation 


Id-) = [l + Go{E+)V] 1$). 


(37) 


This is a convenient approximation that is reliable under 
some limits that include a tight localization in time and 
space of the scatterer, so that we do not need to con¬ 
sider multiple absorptions and reemissions of the scat¬ 
tered particle — a sort of RWA. However, this is not the 
only solution to the problem. If we are able to invert 
the full Hamiltonian approximately we can construct the 
solution |d>) as 

|Mr)=[l + G(F+)I/]|d>) =:!$) +It/.), (38) 

where It/.) is the outgoing or scattered state. To prove 
this expression we use ( p^ and the relations [1 + 
G{z)V][l - Go{z)V] = [1 - Go{z)V][l + G{z)V] = 1. 


B. Single-qubit scattering 

Using the fact that iJ is box-diagonal in the excitation 
number space, we will now compute the resolvent G{E) 


where |d>) was assumed to have a well defined momentum 
ko and thus E'^ = w(fco) + 0'*'. We solve Eq. (35) written 
as, (if+ — E[) lijj) = Grouping terms this leads to 


(U+ - As)cs - gsk-ipk = g*sko (40) 

k 

{E+ - u)k)'ipk - Y 9*skds = 0. (41) 

S 

The second equation may be readily solvedgiving an 
equation for c 

{E+ - As)cs - Y = 9*sko- (42) 

^ ^ -^k 

Note that this equation may be written in a much more 
compact form 


c=(F;+-iL)-Vo, (43) 

where the effective Hamiltonian is the one computed in 
the master equation formalism. For one qubit 

H = A' (44) 


and the scattering process has a resonance at the renor¬ 
malized qubit frequency. For two qubits 


/ A( - 171/2 512 -I- /7i2/2\ 
\gi 2 + * 712/2 A '2 - * 72/2 J 


(45) 


and the scattering resonances appear where predicted by 
the effective two-qubit interaction, broadened by the re¬ 
spective decay rates, 7 and 712 - 

Note that unlike the Markovian equation, which is ap¬ 
proximate, the scattering equations are exact for any sep¬ 
aration of the qubits, because they describe the asymp¬ 
totic states under a stable but sparse flow of individual 
photons. 


C. Relation to transfer matrix models 

The previous subsection shows that the single- and 
many-qubit resonances are renormalized due to the indi¬ 
vidual and collective Lamb shifts. This result seemingly 
contradicts earlier predictions of the scattering mod¬ 
els in one-dimensional waveguides and coupled cavities 
[H [m 113 im, where no such renormalization is evi¬ 
dent. This contradiction is due to the choice of model in 
those earlier works and has to be analized case by case. 

In Refs. |14j and |15j . the coupling is written in position 
space as a boundary condition for two propagating fields. 
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The coupling is local in space at the position of the qubit 
and the spectra of the photon are continued towards uj —>■ 
—oo, introducing two independent propagating fields that 
move left- and right-wards. From the point of view of 
our formalism this results in a constant spectral function 
J{uj). Introducing J{uj) = J in our integrals, together 
with some cut-offs, we recover that the imaginary part 
of the kernel is zero, so that A' = A, and the individual 
qubit renormalization is absent. 

Another popular model is an array of coupled cavities 
modeled by a tight-binding model. This model results 
from applying a rotating-wave approximation in our band 
model, introducing a discrete set of N modes {a„, a\^}n=i 
whose Hamiltonian contains a nearest-neighbor hopping 
term 


N 

= '^ujalan- J ^ (46) 

n—1 {n,m) 


Translational invariance allows us to diagonalize the 
model with a Fourier transform, obtaining 

ujk = u}-2Jcos{k), gsk = (47) 

where now the position Xg is an integer and the coupling 
once more does not depend on the frequency. Integrating 
the kernel for this distribution of frequencies gives 


m = 5 ^ 


\/{J — A w)(t/ + A — cu) 


(48) 


which is real for a qubit inside the band, |A—a;| < J, and 
purely imaginary outside the band. Thus, in the usual 
situation in which the qubit may directly interact with 
the photons. A' = A and we have no Lamb shift. 

The lack of single-qubit and collective Lamb shifts in 
these models is thus a coincidence that arises due to 
the approximations in the photon field and in the qubit- 
photon interaction. It is thus not a physically realistic 
effect and when wishing to do quantitative predictions 
of experiments we should go back to models such as the 
ones in Sect. am 


V. NUMERICAL RESULTS 


We now present two numerical approaches to solve the 
system of two qubits interacting via a transmission line 
in the RWA. The first set of simulations aim at com¬ 
paring with our earlier analytic predictions for the toy- 
model with infinite spectrum ([^. In this case we solve 
the integro-differential equation 0 using the fact that 
we have explicit expressions for the kernel and that we 
have an efficient and stable implicit method to solve the 
equation. 

The second set of simulations are more general in their 
design, working directly with ( 10 ) for a problem with a 


finite size. This is possible because we work with more 
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FIG. 3. Single-qubit spontaneous emission rate, 71 , and fre¬ 
quency renormalization, <5, as a function of the model cut-off, 
ujc in 1^. We plot the numerical result obtained by solving 
(111 with an implicit method (solid), together with the self- 
consistent expression [(|24[ ), dashed] and the usual prediction 
without Lamb shift [(|23[), dotted]. 
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FIG. 4. Two-qubit parameters, 7 , 712 and gi 2 , as a function 
of the model cut-off, ujc in We plot the numerical re¬ 
sult obtained by solving (111 with an implicit method (solid), 
together with our prediction \27\ using the self-consistent so¬ 
lution for A± (blue, dashed) and expressions for 7,712 and 
gi 2 that do not include the Lamb shifts (red, dotted). 


realistic scenarios where the spectrum has a physically 
motivated cut-off, as it is the case of continuous transmis¬ 
sion lines ([^ or photonic crystals ([^. In both cases we 
will solve the single excitation model numerically , con¬ 
firming that many of the features obtained before with 
the analytically solvable models are present. 

The three sections below demand a model fitting mech¬ 
anism that allows us to extract the parameters 7 , 712 , gi 2 
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FIG. 5. Single-qubit spontaneous emission (a), collective decay (b) and interaction (c) for two qubits interacting with a one¬ 
dimensional gapless waveguide, ui{k) = — cos{k6x), as a function of the qubit separation d. As simulation parameters we 

use the same gap for both qubits, A, and the speed of light v as units. The waveguide has a length L = 20Ao and a discretization 
5x = L/N, with varying number of modes: N = 200 (dotted), 400 (dot-dashed), N = 1600 (dashed) and N = 3200 (solid). 


or A' from the dynamics of the qubit, a tool that is also 
essential to analyze the Matrix Product States simula¬ 
tions in Sect. EH 

A. Exponential cut-off 

In this first set of simulations we solve numerically 
equation Our starting point is a formal rewrite 

of the equation 

dtc{t) = — f IC{t — r)e“*^'^c(r)dr, (49) 

Jo 

where we have eliminated free evolution, introducing it 
into the integral. Further integration leaves 

c(T)=c(0)-/' [ /C(r - r)e-*'^^c(r)drdt, (50) 

Jo Jo 

which can also be written 

c{T)=c{0)-[ G{T - t)c{t)dt, (51) 

Jo 

where the function G{t) can be analytically computed 
from ([^, with an expression too large to be reproduced 
here. The last step is to discretize time 

A 

m—1 

+G{tn — tm)c{tm+l)]- (52) 

This equation has been appropriately symmetrized to in¬ 
crease the accuracy and make the formula more stable, 
but it then implies that the equation must be solved im¬ 
plicitly, as c{tn) appears also in the sum. 

Integration formulas apply directly to the problem of a 
single qubit’s spontaneous emission. In Fig. we plot the 
two parameters as obtained from numerical integration 


(solid), and from the Markovian formulas, for a problem 
with a moderate cut-off, uic = lOA and coupling strength 
g = 4%A. Notice that while 5 is exactly approximated by 
the analytical expressions, only the self-consistent renor¬ 
malization approaches the numerical results. 

We can repeat this idea using now two qubits and a 
slightly larger cut-off, Wc = 20A, to ease plotting. The 
results are shown in Fig. Notice how the numeri¬ 
cal simulations (solid) reproduce the oscillations of 7 as 
a function of the distance, which are only captured by 
our self-consistent formulas (dashed). All other schemes, 
such as a partial renormalization ( 22 ) or the usual reso¬ 
nant dipole approximation with no change in A, fail by 
a significant margin (red dashed lines) that worsens with 
the distance. 


B. Finite length transmission line 

We now study the first model for which we do not have 
an analytical solution. This is a discretized model for a 
one-dimensional transmission line or waveguide |19j . The 
waveguide will have a total length L divided into seg¬ 
ments of size Ax. The coupling between these segments 
gives rise to the dispersion relation ([^, which is approx¬ 
imately linear for frequencies smaller than Wc, the cutoff 
frequency 


(jj{k) « ujc\k\5x, (53) 

The speed of light, u = -j^ = Wcdx, will be used to¬ 
gether with the qubit frequency A to adimensionalize all 
quantities: from the unit of length ^ = 1 we obtain the 
wavelength of the qubit Aq = = 27r, and also the 

cutoff frequency 10 ,.= 

For the simulation to actually reproduce propagating 
photons, the waveguide length L also should be larger 
than the qubit wavelength and the total simulation time 
should be shorter than the time for a photon to cross the 
whole line and bounce back (open boundary conditions) 
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or return the qubit from the other side (periodic bound¬ 
ary conditions). In our simulations we take L between 
lOAo and 20Ao and thus T < L/v to avoid the revivals. 

From the dispersion relation and the modes, we solve 
numerically the full equation in the single-excitation limit 
(10). This is done with MATLAB’s adaptive ODE solvers 
using an initial condition in which only one qubit is ex¬ 
cited and the line has no photons. The excited qubit will 
relax and spontaneously emit a single photon which, after 
a finite travel time tfught = d/v through the separation 
d, is scattered by the second qubit. It is important to 
remark the existence of an approximate light cone. Out¬ 
side this light cone, t < tfUght, the first qubit does not 
have enough time to significantly influence its unexcited 
counterpart and thus the Markovian model does not ap¬ 
ply. After the time of flight time tfUght, we may start 
talking about interaction and the collective variables c± 
start having independent dynamics. 


In Fig. 1^ we show all dynamical variables from one 
particular simulation. First of all, notice how |cip -I- 
|c 2 p = |c+p -|- |c_p follows an exponential decay with a 
rate, 7 , that remains approximately constant throughout 
the whole evolution. At t/Hg/it the traveling photon hits 



tA 



FIG. 6 . (a) Qubit variables in a single numerical experiment. 
We plot |c+|^ -I- |c_|^ (solid), |c_p (dashed) and |c+p (dash- 
dot). Note how c+ and c_ depart at the time at which the 
photon arrives to the second qubit. The simulation parame¬ 
ters are L = 207r, N = 800, g — 0.04 and d = Aq. (b) We 
extend the same plot to all other qubit separations, represent¬ 
ing ||c-p — |c+p| as a function of time and distance. Outside 
of the light cone, there is no difference between |c_| and |c+|, 
and the qubits are uncorrelated. As guide to the eye we plot 
the dashed line d = vt. 



FIG. 7. Renormalized gap of the first qubit obtained from the 
approximate period of the functions 312 and 712 , as a function 
of the cut-off ujc, for A = 1, g = 4%A. 


the unexcited qubit, which partially absorbs it. At this 
point, variables c+ and c_ depart and acquire different 
exponential decay rates, whose average is 7 and whose 
differece is 712 . Finally, in Fig. Wp we plot the value of 
|c+p — |c_p as a function of time, thereby showing how 
the light cone extends approximately over all distances 
and how at long times the dynamics is periodic on the 
qubit separation. 

We can simulate multiple different processes to recover 
the dependency of 7 , 712 and gi 2 on the qubit separation 
d. This is shown in figure Fig. for different discretiza¬ 
tions of the chain. In each curve we employed a different 
number of sites N = 200,400,1600,3200, but always the 
same length L = 20Ao. Note that both gi 2 and 712 are 
approximately periodic, as in the analytically solvable 
models, with a period that depends on the discretization 
and thus on the cut-off frequency. This is a sign of the 
qubit gap renormalization. 

Studying the photon wavelength from the zeroes of the 
712 and gi 2 we are able to obtain the dependance of A' as 
a function of 6x or Wc as it is showed in Fig.[^ The value 
of the individual spontaneous emission rate 7 was also 
predicted to oscillate due to the influence of the collective 
Lamb shifts in the decay rate. When we average over 
j{d) over the qubit separation, this collective Lamb shift 
is averaged out and we obtain a mean value that only 
depends on A' and, since J{A') oc A', follows a similar 
curve to Fig. 


C. Photonic crystal 

In this second case instead of having a transmission 
line with free travelling photons, we consider a photonic 
crystal Q. This dispersion relation may be obtained, for 
instance, by coupling together multiple cavities or pat¬ 
terning a waveguide. In this case the photon’s speed, de¬ 
termined by the group velocity, depends on its frequency. 
We will fix ujo = A, so that this group velocity is strictly 
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FIG. 8. Single-qubit spontaneous emission (a), collective decay (b) and interaction (c) for two qubits interacting with a one¬ 
dimensional gphotonic crystal, Lij{k) = wo — Jcos(fc). We use periodic boundary conditions for qubits with gap A = wo, 
N = 800, a waveguide length L = 20Ao and a varying bandwith J/ujo = 0.5 (dotted), 0.6 (dashed), 0.8 (dot-dashed) and 1 
(solid). Note how the period of the curves is not significantly affected by the bandwidth, J. 


the coupling strength J 


Vn = 


doJk 


= J 


(54) 




The actual length of the photonic medium is L = N the 
number of sites in the crystal. However now the dis¬ 
persion relation does not depend on N and the number 
of modes will only influence features such as the prop¬ 
agation time of the photons and the smoothness of the 
curves. 

As before, the numerical integration works with the 
full model ( 10 ), now studying only the influence of the 


band width J and the qubit separation d. Once more, we 
find an approximate light cone, this time governed by the 
group velocity Vg. Once more we compute the collective 
variables c± numerically inside and outside the light cone 
and fit the numerical results to the theoretical curves. 

In Fig. we plot the fitted parameters gi2, 7 and 712 
as a function of the qubit separation for different band- 
widths, J. Note how the spontaneous emission rate and 
the interactions grow with decreasing bandwith. Phys¬ 
ically this may be understood as a consequence of the 
reduced photon velocity and an increase in the spectral 
density around the two-level system, which facilitate the 
interaction between the second qubit and the photon that 
was originally emitted. Note also that due to the very 
small cut-off, which cannot grow beyond 2A, there is a 
negligible renormalization of the qubit frequency, which 
is almost invisible in 7 and impossible to recover from 
the periods of 712 and gi 2 . 


VI. MATRIX PRODUCT STATES 

We finish the theoretical study by performing the same 
simulations in a more general situation, that is the full 
spin-boson model 0 without resorting to the single¬ 
excitation or RWA approximation. We do this by using 
the Matrix Product State ansatz for solving numerically 


the full spin-boson model with both qubits interacting 
with the bosonic bath. The method used is a general¬ 
ization of the one outlined in |19j . introducing one more 
qubit and applying the same model fitting techniques to 
extract the parameters 7 , 712>512 and A'. 

The scope of this paper the study of qubit-qubit inter¬ 
actions for realistic couplings, with g < 5%A, such as the 
ones describing transmon or ordinary flux qubits close to 
an open transmission line. The outcome of these simu¬ 
lations is shown in Fig. where we plot the Markovian 
parameters from one MPS simulation together with a nu¬ 
merical simulation using the RWA and single-excitation 
limit. 

The main message from those simulations is therefore 
that the RWA accurately describes the qubit and pho¬ 
ton dynamics for the weak and strong-coupling regimes. 
The numerical and theoretical results from earlier sec¬ 
tions should therefore apply to ongoing experiments and 
could be verified in state-of-the art setups that already 
study photon-qubit interactions. 


VII. SUMMARY AND CONCLUSIONS 

We have studied the dynamics of one and two qubits 
coupled to an open transmission line. We have related the 
dynamics of these qubits, in a regime of weak or strong 
coupling {g/A < 5%) to the RWA or single-excitation 
equations. From these equations we have derived care¬ 
fully the Markovian limit, studying the approximations 
involved and paying special attention to the single-qubit 
and collective Lamb shifts that arise. We showed how 
the Markovian parameters may carry a significant and 
measurable dependence on the microscopic details of the 
underlying photons, such as their ultraviolet cut-offs and 
the exact shape of the spectral function. Such features 
are present not only on the non-equilibrium dynamics 
of the qubits, but also on the spectroscopy or scattering 
properties. Our predictions have been confirmed for a 
large variety of physical models, and also with numeri- 
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FIG. 9. (a) Comparison between the results of MPS sim¬ 

ulations, 'yl2g^, 712 / 2 ^^ and gw/g^, in solid line, and the 
RWA method, in dashed line, (b) Difference ||c+p — |c_p| 
as a fnnction of time and qubit separation, together with the 
light-cone (dashed). The parameters of the simulation are 
N = 321, g = 0.05, A = 1 and u = 1. 


cal simulations based on Matrix Product States without 
RWA approximations. 

One take-home message of this work is that the mi¬ 
croscopic details of qubits interacting with propagating 
photons can have a measurable impact in the quantum- 
optical properties of the qubits. In order to experimen¬ 
tally probe those effects one needs to be able to change 
such microscopic models. As shown in Ref. |27j . a very 
useful setup to realize such changes has already been 
demonstrated in the lab |28j . The experiment would con¬ 
sist on mobile transmon qubits that are suspended over a 


transmission line, thereby allowing us to change the cou¬ 
pling strength and the qubit-qubit separation, and thus 
obtaining not only the curves A', 7,712 and 312 , but also 
the bare parameters in the absence of renormalizations. 
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project PROMISCE. JJGR acknowledges computer re¬ 
sources and assistance provided by the Centro de Super- 
computacion y Visualizacion de Madrid (CeSViMa) and 
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Appendix A: Master equation constants 


We have related the dynamics of spontaneous emission 
to a linear system of differential equations that is local 
in time. In this setting, the coupling between spins and 
the decay of the excited populations is given by a matrix 

pOO 

Hss'=-i Kss'{T)e''^-'^dT, (Al) 

Jo 

— (Ag A5 i'ys)dss' T iSss* ^ss') 


In order to relate this expression to the scattering equa¬ 
tions we will focus on the time integral, expressing it in 
the original couplings as the limit 




poo 

hm -i / ^ 

^ 0 + Jo . 


= lim V 


k 

9skg*,k 






J 0 




■*5sfc5s'fc 


AC -I- f0+ - Wfc 


(A2) 


Note how this is close to the expression that appears in 
the denominator of the Lippmann-Schwinger equation: 
unlike the master equation, the LS formalism deals with 
stationary states, a situation in which all amplitudes 
and ijjk evolve with the same frequency. In that case we 
have to replace A(, —>■ E, and we recover the effective 
Hamiltonian from the scattering states. 
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